﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using RayDen.Library.Core.Components;

namespace RayDen.Library.Core.Primitives
{
    public interface IWavelengthDependentValue
    {
        float Sample(float lambda);
    }

    public class SpectralDistribution : IWavelengthDependentValue
    {
        protected SpectralDistribution()
        {
            nSamples = 0;
            lambdaMin = lambdaMax = delta = invDelta = 0f;
            samples = null;
        }

        public float[] GetValues()
        {
            return this.samples;
        }

        // samples the SPD by performing a linear interpolation on the data
        public float Sample(float lambda)
        {
            if (nSamples <= 1 || lambda < lambdaMin || lambda > lambdaMax)
                return 0f;

            // interpolate the two closest samples linearly
            float x = (lambda - lambdaMin) * invDelta;
            int b0 = MathLab.Floor2UInt(x);
            int b1 = Math.Min(b0 + 1, nSamples - 1);
            float dx = x - b0;
            return MathLab.Lerp(dx, samples[b0], samples[b1]);
        }

        public void Sample(int n, float[] lambda, float[] p)
        {
            for (int i = 0; i < n; ++i)
            {
                if (nSamples <= 1 || lambda[i] < lambdaMin ||
                    lambda[i] > lambdaMax)
                {
                    p[i] = 0f;
                    continue;
                }

                // interpolate the two closest samples linearly
                float x = (lambda[i] - lambdaMin) * invDelta;
                int b0 = MathLab.Floor2UInt(x);
                int b1 = Math.Min(b0 + 1, nSamples - 1);
                float dx = x - b0;
                p[i] = MathLab.Lerp(dx, samples[b0], samples[b1]);

            }
        }
        public float Y()
        {
            float y = 0f;
            for (var i = 0; i < SampledSpectrum.nCIE; ++i)
                y += Sample(i + SampledSpectrum.CIEstart) * SampledSpectrum.CIE_Y[i];
            return y * 683f;
        }

        public float Filter()
        {
            float y = 0f;
            for (var i = 0; i < nSamples; ++i)
                y += samples[i];
            return y / nSamples;
        }

        public void AllocateSamples(int n)
        {
            samples = new float[n];
        }

        public void FreeSamples()
        {
            samples = null;
        }

        public void Normalize()
        {
            float max = 0f;

            for (var i = 0; i < nSamples; ++i)
                if (samples[i] > max)
                    max = samples[i];

            if (Math.Abs(max) < 0.00001f)
            {
                return;
            }

            float scale = 1f / max;


            for (var i = 0; i < nSamples; ++i)
                samples[i] *= scale;
        }

        public void Clamp()
        {
            for (var i = 0; i < nSamples; ++i)
            {
                if (!(samples[i] > 0f))
                    samples[i] = 0f;
            }
        }

        public void Scale(float s)
        {
            for (var i = 0; i < nSamples; ++i)
                samples[i] *= s;
        }

        public void Whitepoint(float temp)
        {
            var bbvals = new List<double>();

            // Fill bbvals with BB curve
            float w = lambdaMin * 1e-9f;
            for (var i = 0; i < nSamples; ++i)
            {
                // Compute blackbody power for wavelength w and temperature temp
                bbvals.Add(4e-9f * (3.74183e-16f * MathLab.Pow(w, -5f)) / (MathLab.Exp(1.4388e-2f / (w * temp)) - 1f));
                w += 1e-9f * delta;
            }

            // Normalize
            float max = 0f;
            for (var i = 0; i < nSamples; ++i)
                if (bbvals[i] > max)
                    max = (float)bbvals[i];
            float scale = 1f / max;
            // Multiply
            for (var i = 0; i < nSamples; ++i)
                samples[i] *= (float)(bbvals[i] * scale);

        }

        public RgbSpectrum ToRgb()
        {
            var c = new RgbSpectrum();
            float yInt = 0f;
            for (var i = 0; i < SampledSpectrum.nCIE; ++i)
            {
                float v = Sample(i + SampledSpectrum.CIEstart);
                yInt += SampledSpectrum.CIE_Y[i];
                c.c1 += v * SampledSpectrum.CIE_X[i];
                c.c2 += v * SampledSpectrum.CIE_Y[i];
                c.c3 += v * SampledSpectrum.CIE_Z[i];
            }
#if !VERBOSE
            Assert.IsNotNaN(c.c1);
            Assert.IsNotNaN(c.c2);
            Assert.IsNotNaN(c.c3);
#endif

            var xyz = (c / 683f);

            float[] rgb;
            SampledSpectrum.XYZToRGB(xyz.ToArray(), out rgb);
            c = new RgbSpectrum(rgb);
            //res.GammaCorrection();
            return c;
        }


        public RgbSpectrum ToXyz()
        {
            var c = new RgbSpectrum();
            for (var i = 0; i < SampledSpectrum.nCIE; ++i)
            {
                float v = Sample(i + SampledSpectrum.CIEstart);
                c.c1 += v * SampledSpectrum.CIE_X[i];
                c.c2 += v * SampledSpectrum.CIE_Y[i];
                c.c3 += v * SampledSpectrum.CIE_Z[i];
            }

            var xyz = (c / 683f) * SampledSpectrum.yint;
            return xyz;
        }

        protected int nSamples;
        protected float lambdaMin, lambdaMax;
        protected float delta, invDelta;
        protected float[] samples;
    }

    public class RegularSPD : SpectralDistribution
    {
        public RegularSPD() : base() { }
        public RegularSPD(float[] s, float lMin, float lMax, int n)
        {
            init(lMin, lMax, s, n);
        }
        protected void init(float lMin, float lMax, float[] s, int n)
        {
            lambdaMin = lMin;
            lambdaMax = lMax;
            delta = (lambdaMax - lambdaMin) / (n - 1);
            invDelta = 1f / delta;
            nSamples = n;

            AllocateSamples(n);

            // Copy samples
            for (var i = 0; i < n; ++i)
                samples[i] = s[i];
        }
    }

    public enum SPDResamplingMethod { Linear, Spline };

    public class IrregularSPD : SpectralDistribution
    {
        public IrregularSPD() : base() { }
        public IrregularSPD(float[] wavelengths, float[] samples,
                            int n, float resolution = 5,
                            SPDResamplingMethod resamplignMethod = SPDResamplingMethod.Linear)
            : base()
        {


            float lambdaMin = wavelengths[0];
            float lambdaMax = wavelengths[n - 1];

            var sn = MathLab.Ceil2UInt((lambdaMax - lambdaMin) / resolution) + 1;

            List<float> sam = Enumerable.Range(0, sn).Select(item => 0f).ToList();


            if (resamplignMethod == SPDResamplingMethod.Linear)
            {
                var k = 0;
                for (var i = 0; i < sn; i++)
                {
                    float lambda = lambdaMin + i * resolution;

                    if (lambda < wavelengths[0] || lambda > wavelengths[n - 1])
                    {
                        sam[i] = 0f;
                        continue;
                    }

                    for (; k < n; ++k)
                    {
                        if (wavelengths[k] >= lambda)
                            break;
                    }

                    if (wavelengths[k].NearEqual(lambda))
                        sam[i] = samples[k];
                    else
                    {
                        float intervalWidth = wavelengths[k] - wavelengths[k - 1];
                        float u = (lambda - wavelengths[k - 1]) / intervalWidth;
                        sam[i] = MathLab.Lerp(u, samples[k - 1], samples[k]);
                    }
                }
            }
            else
            {
                List<float> sd = new List<float>(n);

                calc_spline_data(wavelengths, samples, n, sd);

                var k = 0;
                for (var i = 0; i < sn; i++)
                {
                    float lambda = lambdaMin + i * resolution;

                    if (lambda < wavelengths[0] || lambda > wavelengths[n - 1])
                    {
                        sam[i] = 0f;
                        continue;
                    }

                    while (lambda > wavelengths[k + 1])
                        k++;

                    float h = wavelengths[k + 1] - wavelengths[k];
                    float a = (wavelengths[k + 1] - lambda) / h;
                    float b = (lambda - wavelengths[k]) / h;

                    sam[i] = Math.Max(a * samples[k] + b * samples[k + 1] +
                        ((a * a * a - a) * sd[k] + (b * b * b - b) * sd[k + 1]) * (h * h) / 6f, 0f);
                }
            }

            init(lambdaMin, lambdaMax, sam, sn);
        }

        protected void init(float lMin, float lMax, IList<float> s, int n)
        {
            lambdaMin = lMin;
            lambdaMax = lMax;
            delta = (lambdaMax - lambdaMin) / (n - 1);
            invDelta = 1f / delta;
            nSamples = n;

            AllocateSamples(n);

            // Copy samples
            for (var i = 0; i < n; ++i)
                samples[i] = s[i];
        }

        protected void calc_spline_data(float[] wavelengths, float[] amplitudes, int n, IList<float> spline_data)
        {
            List<float> u = new List<float>(n - 1);

            // natural spline
            spline_data[0] = u[0] = 0f;

            for (var i = 1; i <= n - 2; i++)
            {
                float sig = (wavelengths[i] - wavelengths[i - 1]) /
                    (wavelengths[i + 1] - wavelengths[i - 1]);
                float p = sig * spline_data[i - 1] + 2f;
                spline_data[i] = (sig - 1f) / p;
                u[i] = (amplitudes[i + 1] - amplitudes[i]) /
                    (wavelengths[i + 1] - wavelengths[i]) -
                    (amplitudes[i] - amplitudes[i - 1]) /
                    (wavelengths[i] - wavelengths[i - 1]);
                u[i] = (6f * u[i] /
                    (wavelengths[i + 1] - wavelengths[i - 1]) -
                    sig * u[i - 1]) / p;
            }

            float qn, un;

            // natural spline
            qn = un = 0f;
            spline_data[n - 1] = (un - qn * u[n - 2]) /
                (qn * spline_data[n - 2] + 1f);
            for (var k = 0; k < n - 1; ++k)
                spline_data[n - k - 1] = spline_data[n - k - 1] *
                    spline_data[n - k] + u[k];
        }

    }
}
